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We have studied the time evolution of the heavy ion luminosity and bunch intensities in the 
Relativistic Heavy Ion Collider (RHIC), at BNL, and in the Large Hadron Collider (LHC), at 
CERN. First, we present measurements from a large number of RHIC stores (from Run 7), collid- 
ing 100 GeV/nucleon 197 Au 79+ beams without stochastic cooling. These are compared with two 
different calculation methods. The first is a simulation based on multi-particle tracking taking into 
account collisions, intrabeam scattering, radiation damping, and synchrotron and betatron motion. 
In the second, faster, method, a system of ordinary differential equations with terms describing 
the corresponding effects on emittances and bunch populations is solved numerically. Results of 
the tracking method agree very well with the RHIC data. With the faster method, significant dis- 
crepancies are found since the losses of particles diffusing out of the RF bucket due to intrabeam 
scattering are not modeled accurately enough. Finally, we use both methods to make predictions of 
the time evolution of the future 208 p D S2 + beams in the LHC at injection and collision energy. For 
this machine, the two methods agree well. 

PACS numbers: 29.20.db, 25.75.-q 



I. INTRODUCTION 

During the design and operation of a heavy-ion col- 
lider, e.g., the Relativistic Heavy Ion Collider (RHIC) |l| 
or the Large Hadron Collider (LHC) [2|, one of the main 
goals is to maximize the time-integral of the luminosity 
C at the experiments. As usual, luminosity is defined as 
the event rate per unit cross section and depends only on 
the spatial density distributions of the colliding beams at 
the collision points. 

The time-evolution of the spatial densities is deter- 
mined via single-particle dynamics from that of the 
phase-space distributions, i.e., the kinetics of the beams. 
This, in turn, is determined by the combined influences 
of several inter-dependent physical processes. The ac- 
tion of some of these (e.g., beam-beam collisions, scat- 
tering on residual gas) principally remove particles from 
the beams. Others (e.g., intrabeam scattering (IBS), ra- 
diation damping) predominantly change their distribu- 
tion in space and momentum. To maximize J C dt, losses 
caused by other processes than the collisions themselves 
need to be minimized and the beam sizes should stay 
small. These processes therefore need to be understood 
and modeled in quantitative detail. 

There have been numerous previous studies of the time 
evolution of luminosity (in colliders) or other perfor- 
mance parameters such as the emittances (e.g., in syn- 
chrotron light sources or the damping rings of linear col- 
liders). Among many possible references, we cite 



which are particularly close to the applications of this 
paper. 

Many of these are based on the solution of systems 
of coupled ordinary differential equations (ODEs) that 
describe the time evolution of a few parameters charac- 
terizing the beam distributions, typically the intensities 
and the first- and second-order moments of the distribu- 
tions (beam centroids and emittances) . Such systems can 
only be closed with additional assumptions on the nature 
of the beam distributions. Typically these are assumed 
to be gaussian in all three degrees of freedom. Besides 
closing the system of equations, this can also allow con- 
venient analytical forms for some of the terms, e.g., those 
describing intrabeam scattering or the way in which the 
luminosity is modified by crossing angles at the interac- 
tion point. This approach was applied to the evolution of 
heavy- ion luminosity in the LHC in Ref. @ . Alternative 
approaches involve solving the Fokker-Planck equation 
for the beam distribution functions 0, Q and single par- 
ticle tracking. In Ref. j| both the ODE method and 
particle tracking were implemented in the same simula- 
tion code. 

In this article, we study the time evolution of colliding 
heavy-ion beams in RHIC and LHC. In Sec. In] we present 
measurements of luminosity and beam intensities during 
a large number of physics stores during Run-7 
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RHIC and give a summary of the running conditions. 
The data are compared with simulations based on two 
different models: 

• Multi-particle tracking: This method, described in 
Sec. IIII1 is rather direct and accurate but slow. We 
use an extended version of the code introduced in 
Refs. motivated by the fact that the longi- 
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TABLE I. Typical beam and machine parameters for RHIC Run-7 and the LHC, given for the beginning of store. The transverse 
emittance in RHIC showed large variations between stores [Tc( |. 



Parameter 


Unit 


RHIC collision 


LHC collision 


LHC injection 


Species 




lav a va+ 
Au 






Beam energy 


GeV /nucleon 


100 


2759 


177.4 


Lorentz factor 7 rc i 




107.4 


2963.5 


190.5 


Bunch intensity Ni, 


10 H 


1.1 


0.07 


0.07 


Bunches per beam 




103 


592 


592 


Normalized transverse rms emittance 


fim 


3.1 


1.5 


1.4 


Long, rms emittance 


eV s/nucleon 


0.25 


0.25 


0.07 


rms bunch length 


cm 


30 


7.94 


9.97 


i uin enei gy api eatr 


1 n _a 


n s 

u.o 


nil 


u.oy 


No. active interaction points (IPs) 




2-4 


1-3 





Optical function at IP f$* xv 


m 


0.8 


0.5-0.55 




Crossing angle at IP 


/^rad 





70-285 




Peak luminosity 


10 27 cm^s- 1 


3 


1 




RF harmonic numbers h 




360, 2520 


35640 


35640 


RF gap voltage 


MV 


3 (h = 2520) 
0.3 (ft = 360) 


16 


8 



tudinal bunch distributions in RHIC are not gaus- 
sian, as discussed in Sec. IIIIEI It is important to 
have a simulation that can model IBS for arbitrary 
profiles of bunched beams. 

• ODEs: Numerical solution of a coupled system of 
ODEs describing the beam emittances and bunch 
populations as in Ref. @. This method, discussed 
in Sec. lIVi is fast but lacks accuracy when the beam 
distributions are not gaussian. 

In Sec. |V[ both methods are compared with the data 
from RHIC, and in Sec. I VII we make predictions for the 
LHC. 

II. MEASURED LIFETIMES IN RHIC 

The RHIC Run-7 consisted of 191 physics stores in 
2007 with colliding 197 Au 79+ beams at an energy of 
100 GeV/nucleon. The main beam and machine param- 
eters for Run-7 are summarized in Table HI 

During Run-7, longitudinal stochastic cooling became 
operational for the Yellow ring [14] (the two rings of 
RHIC are called Blue and Yellow). Stochastic cooling 
counteracts the longitudinal diffusion caused by intra- 
beam scattering (IBS) that eventually pushes particles 
outside the longitudinal acceptance (RF bucket). This 
loss process, discussed in Sec. MI El is responsible for a 
large fraction of the beam losses in RHIC. In this article, 
we focus on stores without stochastic cooling in order to 
make a comparison with the LHC. 

The time-dependent average bunch intensities N% (i = 
1, 2 for Blue and Yellow respectively) were measured us- 
ing DC transformers [l[ and can be fitted empirically by 
an exponential function 

N m (t) = N(0)exp(-t/T N ). (1) 



The beam lifetimes T/v fitted to data from the first 3 h of 
all physics stores are shown in Fig. Q] and the mean and 
standard deviation values are listed in Table [TTJ Reg- 
ular physics stores last 5 h, but some stores terminate 
prematurely. Almost all stores have data for 3 h. 

The luminosity was recorded during every store by 
detecting neutrons emitted by ions that have under- 
gone mutual electromagnetic dissociation in the colli- 
sions [H, [n| . It turns out that it can be well fitted by a 
sum of fast- and slow-decaying exponential functions 

£ fit (t) = Aexp(-t/T f ) + Bexp(-t/T s ) (2) 

with the fit parameters A,Tf,B,T s . This is a 3- 
parameter fit since £gt(0) = A + B. The fit function 
@ is purely phenomenological: it is not chosen on the 
basis of any particular physical model but rather for the 
fact that it generally fits rather well and is convenient to 
implement. The fit parameters for all stores are shown in 
Fig. [5] and the means and standard deviations in Table [III 
An example of the measured and fitted luminosity and 
bunch intensity is shown in Fig. [3] 

III. PARTICLE TRACKING SIMULATION 

The tracking program is based on a code initially used 
for stochastic cooling [l2T[T^ |. which has been developed 
further to simulate the time evolution of two circulating 
and colliding bunches. 

Each bunch contains a number of macro particles (in 
the simulations, 5x 10 4 were used), which are tracked in a 
6D phase space. Normalized coordinates (x,p x ,y,p y ) are 
used in the transverse planes and (t,pt) with p t — 7 — 70 
in the longitudinal (70 is the Lorentz factor of the refer- 
ence particle) . The particle coordinates are updated on a 
turn-by-turn basis as a function of the physical processes 
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TABLE II. Average and rms values of fitted lifetime param- 
eters for all 191 physics stores of RHIC Run-7. 



Parameter 


Unit 


Average 


Standard 
deviation 


Blue beam lifetime Tjv 


h 


9.9 


1.7 


Yellow beam lifetime Tjv 


h 


10.9 


3.3 


Luminosity, fast decaying 
component A/(A + B) 


% 


33 


12 


Luminosity lifetime T/, 
fast decaying 


h 


0.6 


0.5 


Luminosity, slow decaying 
component B/(A + B) 


% 


67 


12 


Luminosity lifetime T s , 
slow decaying 


h 


3.9 


2.5 



Blue 



T N (h) 
Yellow 






A/(A+B) 




0.6 0.! 

B/(A+B) 




7}(h) 



2 3 4 

r s (h) 



FIG. 1. Blue and Yellow beam lifetimes Tm obtained by fit- 
ting Eq. |T} to measurements. Only stores without stochastic 
cooling are included. 



FIG. 2. Distributions of fitted luminosity lifetime parameters 
defined by Eq. where T/ is chosen to be the faster decay 
time and T 3 the slower. 



acting on them. The following processes are taken into 
account: 

• Synchrotron and betatron motion: All coordinates 
are updated deterministically on every turn, taking 
into account the machine tune, chromaticity and 
RF voltage. Particles outside the time acceptance 
of the RF bucket are removed. In RHIC and LHC, 
these coasting particles are cleaned by extracting 
them resonantly through a few dipole kicks, send- 
ing them onto the collimators. Therefore they have 
a negligible influence on beam dynamics. Linear 
betatron coupling is taken into account in thin- 
lens approximation. A detailed description can be 
found in Ref. [Hj]. 

• Collisions: collision probabilities are sampled as a 
function of the local density of the opposing beam 
at (x, y). This is described in Sec. IIII A[ 

• IBS: Each particle is given a random kick, calcu- 
lated using the Piwinski model, but modulated by 



the local particle density in order to account for 
non-gaussian longitudinal bunch profiles. This is 
described in Sec. IIII CI 

• Radiation damping and quantum excitation: All 
particles receive a deterministic amplitude decay 
and a random excitation. This is described in 
Sec.IUDj 

Other processes are neglected. Beam loss rates in RHIC 
do not change visibly when the beams are brought into 
collision, and the beam-beam parameter is even smaller 
for ions in the LHC [2j. Therefore, the beam-beam ef- 
fect is neglected in our simulations. A detailed treat- 
ment of beam-gas scattering has been done elsewhere for 
RHIC 0,U| and LHC In both cases, lifetimes 

of the order of 100 or several hundreds of hours were 
found. The emittance blowup at RHIC is fractions of a 
percent per hour [l8[ and standard formulas [2l[ give a 
similar result for the LHC with predicted vacuum condi- 
tions P. Il9j]. Therefore, beam-gas scattering has a neg- 
ligible impact on the dynamics in RHIC and LHC when 
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FIG. 3. A typical example of measured and fitted luminosity 
(top) and bunch population (bottom). The fits are done only 
over the first three hours, indicated by the dashed vertical 
line, as some stores terminate prematurely. The shown fits 
have x 2 /DOF = 1.2 and 0.86. The measurement errors were 
assumed to be 1.4 x 10 25 cm~ 2 s _1 and 2.5 x 10 6 respectively, 
given by the spread of measurement points over a short time 
interval. 



compared to collisions or IBS. 

Furthermore, calculations with the mad-x [22] pro- 
gram show that the Touschek effect (large-angle scat- 
tering from the center of the bucket bringing particles 
outside the energy acceptance) is negligible in both ma- 
chines, with lifetimes of hundreds or thousands of hours. 
However, small-angle scattering of particles already close 
to the acceptance is important and included as discussed 
in Sec. lilTO 

The strengths of the processes (collision probabilities, 
inverse radiation damping times, and kicks from intra- 
beam scattering and quantum excitation) are scaled up 
to account for both the smaller number of macro particles 
and by an additional factor, set by the user, that reduces 
the computational time, so that the number of turns in 
the simulation corresponds to a much larger number of 
turns in the real machine. Typically one simulation turn 
corresponds to 2 x 10 4 real machine turns in our simu- 
lations to achieve a relatively fast execution without loss 
in precision. 



A. Collisions 

A collision probability P\ is calculated for every par- 
ticle on every turn, and a random number is sampled to 
determine if an interaction takes place. In that case, the 
particle is removed. To calculate P\ we perform an over- 



lap integral of the density of the opposing bunch with a 
Dirac (5-function that represents the trajectory of a single 
particle. The details of this are shown in Appendix [X] 
The most general collision routine makes no assumptions 
on the beam distributions and uses a discrete binning of 
the particles. The integration is then replaced by a sum 
over the bins. This routine is slow and a much faster 
code is obtained with some simplifications. 

First we assume gaussian distributions in x and y, 
which is generally a very good approximation (in RHIC, 
the longitudinal distribution is however non-gaussian as 
explained in Sec. IIII E|) . Furthermore, we assume a com- 
mon luminosity reduction factor instead of computing it 
for every particle, so that P\ is a function of the trans- 
verse coordinates only. In the transverse plane we use 
action-angle variables {J X i4>x) defined for a particle in 
bunch 1 by 



Xi = ^2J x /3* y cos (f> x 



l2Jx 

P*xy 



a xy cos 



<t>x) 1 



(3) 



with the angle variables given at the IP and an analo- 
gous definition in y. Here f3* y is the value of the optical 
function at the IP (assumed to be the same in x and y) 
where we define the longitudinal coordinate to be s = 
and x' = dx/ds. We note that a* xy = -f3' xy (0)/2 = 0. 

Since every simulation turn corresponds to a large 
number of machine turns, where cj> has a uniform distri- 
bution on the interval [0, 2n], we average Pi over (j>. The 
resulting collision probability for a particle in bunch 1, 
derived in Appendix El is a function of the betatron ac- 
tions J only: 



Pi = aN' 



Ri, 



(4) 

Here a is the interaction cross section, N 2 the intensity 
of bunch 2, Iq a modified Bessel function, (e.2x,£2y) are 
the geometric rms emittances of bunch 2 and i?tot is the 
global luminosity reduction factor, which takes into ac- 
count the hourglass effect and a crossing angle 28. Writ- 
ing the longitudinal distributions of the two bunches, nor- 
malized to unity, as p Z i{zi), with being the distance to 
the center of bunch i and i — 1,2, the reduction factor is 
given by 



exp - 



2/3* y sm 2 



2/3^ 2 [l+cos(2e)] 



(zi+z 2 ) 
Pzl(zi)p z2 (z2) 



1 



dzidz2. (5) 



s 2 e 



When collisions occur at several IPs, Pi is summed up 
over all of them (only the factor R tot / j3* y varies). Fur- 
thermore, Pi has to be scaled to account for the smaller 
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number of macro particles and the shorter time scale in 
the simulation. 

If the transverse action in bunch 1 is assumed to be dis- 
tributed as p\u — exp(— Ju/ei«)/ei« for u = x,y, equiva- 
lent to a gaussian distribution in u\ and p u \, we calculate 
the luminosity as 



£ = k h f TCV Ni / pi x piy 



fcb/rev^l^2 



i dJ x dJ 11 



X 



i?t, 



(6) 



where fcb is the number of bunches circulating at fre- 
quency / r ev in each ring. Eq. (|6]) agrees with standard 
formulas I2ll. 



B. Interaction cross section 

To calculate P\ we need the total cross section a for 
interactions in which colliding particles are lost from the 
beam. Apart from inelastic hadronic interactions, bound- 
free pair production (BFPP) and electromagnetic disso- 
ciation (EMD) have to be taken into account. These 
electromagnetic processes create particles with a charge- 
to-mass ratio different from the main beam, which even- 
tually leads to particle loss. Detailed discussions of these 
loss mechanisms in heavy ion colliders can be found in 
Refs. 

Cross sections were calculated for BFPP in Ref. [26[ 
and for EMD in Ref. [27j . We have taken standard values 
for the inelastic cross sections as used by the RHIC and 
LHC experiments @, I281 - I301 ] . In Table IIIII we summarize 
the cross sections of the different processes that were used 
in the simulation. As can be seen, the majority of the 
luminosity is used up by BFPP and EMD rather than 
the inelastic interactions, which are the usual object of 
study of the experiments. 



TABLE III. Interaction cross sections a in RHIC and LHC 
for different processes removin g io ns from the beam. Values 
are taken from Refs. @, [M E3,lll • The EMD cross sections 
include all decay channels. 



Process 


a in RHIC (barn) 


a in LHC (barn) 


BFPP 


117 


281 


EMD 


99 


226 


Inelastic 


7 


8 


Total 


223 


515 



C. Intrabeam scattering 

The standard formalisms for describing IBS [3ll-[36l| as- 
sume gaussian profiles. To simulate RHIC, where the 
longitudinal distributions are known to be non-gaussian, 



an IBS model has to go beyond this assumption. Such a 
model has been implemented in Refs. [H, HH, which we 
summarize here for completeness. It assumes that the 
longitudinal and transverse degrees of freedom are in- 
dependent and that the transverse distributions remain 
gaussian. 

For simplicity and speed, the IBS routine starts from 
the Piwinski model [3l[. Our code gives the user the 
choice between a smooth lattice approximation, result- 
ing in fast execution, or a slower but more precise calcu- 
lation taking the full lattice into account. In this case, 
the optical functions are read in from an external ta- 
ble created by MAD-x. Furthermore, the IBS rise times 
can be calculated with the term n 2 /(3 u (i] being the dis- 
persion and fi u the optical lattice function) replaced by 
% = [rj 2 + (j3 u r]' — |/3„J?)] as proposed in Ref. [35[. Here 
primed quantities represent derivatives with respect to 
the path length s. We call this model the modified 
Piwinski method and the model presented in Ref. |3l| 
the original Piwinski method. We define the rise times 
Tms,u,u = x,y,z by 



di€u 

dt Tibs,u 



(7) 



The rise times are calculated for gaussian distributions 
but are then modulated for every particle by the local 
beam density to account for arbitrary longitudinal pro- 
files. Thus, particles are given normally distributed mo- 
mentum kicks A.p u in each plane on every simulation 
turn: 



Ap u = ra m J2T 1B 1 s ^T rev a t y/Trpt(t) 



(8) 



Here r is a gaussian random number with zero mean and 
unit standard deviation, T rcv the revolution time, at the 
standard deviation of t and a pu the standard deviation 
of the momentum p u in plane u. It can be shown that if 
the longitudinal distribution p t , normalized to unity, is 
gaussian, integrating over all particles yields exactly the 
growth given by Eq. (JT)). The strength of the IBS kicks 
are also scaled up by the time ratio of the simulation. 

The rise times at the beginning of a store, calculated 
with the different IBS models in the tracking code and 
parameters in Table HI are shown in Table IIVI We show 
also results from MAD-x, as discussed in Sec. IIVI and, for 
comparison, the rise times calculated using the Bjorken- 
Mtingwa method [32| (with the Coulomb logarithm cal- 
culated by mad-x). As can be seen, IBS is strong in 
RHIC and at injection energy in the LHC. 

For the cases studied here, the Piwinski models tak- 
ing the lattice into account agree well with the Bjorken- 
Mtingwa method, which is considered more general [35] . 
The smooth lattice approximation, with the average j3- 
function calculated from the tune, gives significantly 
shorter rise times, thus overestimating the strength of 
IBS. The average difference between the modified Pi- 
winski method, taking the lattice into account, and the 
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Bjorken-Mtingwa method is 2.9% and the maximum dif- 
ference 5.1%. Despite some known deficiencies in mad- 
X [37j, the code agrees well with the other methods for 
the cases studied here. 

Because of the better agreement with the other more 
general methods, we use the modified Piwinski method 
for all simulations unless otherwise stated. We take the 
found discrepancies as a guide to the maximum uncer- 
tainty of the rise times used in the simulation. The in- 
fluence of the uncertainty of the IBS model on the time 
evolution of the luminosity and bunch population is very 
small compared to other error sources. A comparison 
between two of the models is shown in Appendix [Bl 

TABLE IV. IBS rise times for the emittances in RHIC and 
LHC, calculated using the original Piwinski model in smooth 
lattice approximation (Piw. sm.), original Piwinski taking the 
lattice into account (Piw. latt.), and the modified Piwinski 
method including the lattice (mod. Piw.). The calculations 
for RHIC are done for gaussian bunches — these values are 
modulated for each particle by the local density in the track- 
ing code. Results from the Bjorken-Mtingwa (B-M) method 
and MAD-x, used in Sec. IIVI are also shown. 





RHIC coll. 


LHC coll. 


LHC injection 


Tms,x, mod. Piw. (h) 


1.86 


13.8 


6.15 


Tibs.x, Piw. latt. (h) 


2.05 


14.8 


6.62 


Tms,x, Piw. sm. (h) 


1.67 


10.3 


4.42 


Iras,*, MAD-x(h) 


2.06 


13.2 


6.77 


Tibs,*, B-M (h) 


1.88 


14.1 


5.98 


Tibs.z, mod. Piw. (h) 


1.96 


8.26 


2.84 


Tibs,*, Piw. latt. (h) 


1.93 


8.14 


2.82 


Tibs,*, Piw. sm. (h) 


1.50 


7.52 


2.50 


Tibs,*, MAD-x(h) 


2.03 


7.89 


3.03 


Tibs,*, B-M (h) 


1.99 


7.84 


2.78 



D. Radiation damping 

Radiation damping and quantum excitation are mod- 
eled in each plane u by a deterministic decay, given by the 
emittance damping time T ra d iU , and a random excitation, 
as described in Ref. 38]. Since T radiU ^> T rcv , we expand 
to first order in T re v/T ral ± tU , so that the decay coefficient 
for one turn becomes exp(-T rcv /Trad, n) ~ l-T r c V /T radj „. 
The quantum excitation is determined by the usual radi- 
ation integrals over the lattice [ll[ and would lead, in the 
absence of other dissipative effects, to stationary gaussian 
distributions with rms sizes a c(i . u . The corresponding 
one-turn map for the momentum coordinate p t is there- 
fore 

p t -> p t (l - T Tev /T TaidjZ ) + cr ecitZ rJ 2(r 2 )T rov /T radiZ . (9) 

where r is a unit, zero-mean random number. Similar 
maps are used in the transverse planes. 



In both RHIC and LHC, the photons emitted by heavy 
nuclei are very soft |2| so the cr q,u are several orders of 
magnitudes smaller than the real beam and quantum ex- 
citation has no significant effect on the dynamics. In 
RHIC, radiation damping is also negligible but it is im- 
portant, enough to counter IBS, in the LHC at collision 
energy 0, Q . The damping times are summarized in Ta- 
bleEl 



TABLE V. Calculated radiation damping times for the emit- 
tances in RHIC and LHC. The damping is equal in the trans- 
verse planes. 





RHIC collision 


LHC collision 


LHC injection 


Trad, z (h) 


825 


6.3 


23749 


T Ta ,d,xy (h) 


1650 


12.6 


47498 



E. Starting distribution 

It is important to choose an appropriate, realistic, 
starting distribution of the particles in the tracking; oth- 
erwise agreement between measurement and simulation is 
poor. In the transverse plane, the beams are well approx- 
imated by gaussian distributions, both in RHIC and (it 
is expected) LHC, and the initial coordinates are gener- 
ated by the code from the starting emittances. The data 
from the RHIC stores does not in fact include the mea- 
sured transverse beam sizes so these have been inferred 
from the logged initial bunch populations and luminosity 
using Eq. ([6]). A strong coupling between the horizontal 
and vertical planes, keeping the beams round on a short 
time scale, was assumed. This has been observed em- 
pirically in RHIC fl2j , while for the future ion operation 
in the LHC, this corresponds to the somewhat idealized 
situation described by Ref. [3^. We also assumed equal 
transverse beam sizes in every bunch of the two beams. 
Despite not being justified by measurements, this reduces 
the dimensionality of the problem to a tractable level; 
there is, in any case, insufficient data to attempt a fit 
to the small variations in intensity and size of individual 
bunches. Indeed we have found that, if different beam 
sizes are used in the simulation, the luminosity remains 
approximately unchanged while the bunch currents show 
small variations. 

In the LHC, the bunch length is much shorter than 
the RF bucket size and an approximation of a gaussian 
distribution was used. This is not the case in RHIC, 
as can be seen in Fig. 21 where an example of the bunch 
current measured at RHIC when the beams were put into 
collision is shown. To understand the bunch shape, we 
discuss briefly the synchrotron motion in RHIC. 

RHIC uses a double RF system (see Table H]) and the 
longitudinal emittance is comparable to the bucket size 
of the h = 2520 system. The combined RF voltage from 
both systems is shown in Fig. [5] We show also its neg- 
ative time integral, which is proportional to the "poten- 
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-15 -10 -5 5 10 15 

t(ns) 

FIG. 4. Example of a measured bunch profile in RHIC after 
the ramp when collisions are turned on. We show also the 
simulated profile used as starting conditions in the simula- 
tions. It was generated by letting an initially gaussian bunch 
circle in the machine under the influence of IBS. 




-40 -30 -20 -10 10 20 30 40 

t(ns) 




-40 -30 -20 -10 10 20 30 40 
t(ns) 



FIG. 5. The RF voltage (parameters given in Tabl^J) and its 
negative time integral. Since synchrotron radiation is negligi- 
ble, the synchronous phase is at the origin at Vrf ~ 0. 



system is turned on. To make the bunches fit in the 
h = 2520 buckets, they are shortened through a rotation 
in the longitudinal phase space, but it is inevitable that 
some particles escape into neighboring buckets [4Cj |. 

To reproduce the profile in Fig. @] and keep the phase 
space consistent with the RF motion, we used as starting 
conditions coordinates generated by tracking an initially 
gaussian bunch, located only in the central h = 2520 
bucket, under the influence of IBS only. In Fig. [5] we 
show the longitudinal phase space from this simulation 
on different turns. Only particles with an amplitude high 
enough to pass from the central h = 2520 bucket are 
present in the neighboring buckets. Therefore, the cen- 
ters of these buckets are empty in phase space. This gives 
a time profile similar to the expected result of the bunch 
shortening (see Ref. [40l|). 

Fig. 2] shows the resulting bunch current for the popu- 
lation that was used as starting distribution in the RHIC 
simulations. It corresponds well to the measured profile 
except around t = ±5 ns, at the centers of the first side 
buckets (see Fig. [SJ- The discrepancy means that in the 
measurement, the phase space distribution is less hollow 
than in Fig. [5] 

Even though it would be possible to recreate the mea- 
sured time profile exactly, this ad-hoc construction would 
depend on an arbitrary choice of the energy deviations of 
the extra particles in the side buckets. We prefer to keep 
the starting conditions as simple as possible in order to 
test the predictive power of the code for a large number 
of different conditions and therefore use the profile shown 
in Fig. [U 

In Fig. [7] we show an example of the simulated losses 
caused by IBS and collisions in RHIC. Simulations show 
that if no particles are present in the outer h = 2520 
buckets in the starting population, losses from debunch- 
ing do not occur from the beginning of the store. With 
such a longitudinal starting distribution, the initial slope 
of Ni(t) is less steep and the agreement with measure- 
ments is significantly worse. 

It should be noted that the assumptions discussed in 
this section fully determine all free parameters in the 
simulation. 



tial energy" term in the Hamiltonian of the synchrotron 
motion. One can thus picture the particles "sliding" on 
this surface. As a particle in the central bucket con- 
tinuously receives small momentum kicks from IBS, it 
oscillates with higher amplitudes until it has enough en- 
ergy to enter the next h — 2520 bucket. When it has 
a high enough amplitude to leave the h = 360 bucket, 
defining the acceptance, the motion becomes unbounded 
and it is considered lost by the simulation. This loss 
process is called debunching and is very strong in RHIC. 
For this reason, very significant improvements of the life- 
times and integrated luminosity were possible through 
the implementation of stochastic cooling fl2l - [l4j |. 

The profile in Fig. @] results from the RF gymnastics 
performed after the energy ramp, when the storage RF 



IV. ODE SIMULATION 

As an alternative to the relatively slow tracking sim- 
ulation we compare it with a faster method, which does 
not follow single particles but the RMS emittances. This 
is done by solvingnumerically a system of coupled ODEs, 
similar to Refs. ](| [t| |4l|. In all cases treated here, we 
assume no coherent oscillations of the bunches so that 
the first-order moments vanish. Then we assume both 
beams to be gaussian in all degrees of freedom, which 
obviously is not true for the longitudinal plane in RHIC, 
but is expected to hold in the LHC. Solutions for RHIC 
are presented in Sec. |V] anyway for comparison. 

We assume round beams (e X i = e y i = e xy i for beams 
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t = min 




t(ns) 
t = 9 min 




t(ns) 
t = 67 min 




t(ns) 



FIG. 6. A simulated bunch under influence of IBS only in 
longitudinal phase space at different times. Each point rep- 
resents one particle out of 50000. The center of the bunch 
is located in the central bucket at t = 0. The longitudinal 
momentum is defined as pt = 7 — 70 , where 70 is the Lorentz 
factor of the synchronous particle. The figure shows an ideal- 
ized situation, where all particles start in the central bucket 
(reality is different because of the RF gymnastics). The dis- 
tribution at 67 minutes was used as starting conditions in the 
longitudinal plane for the full tracking, due to its similarity 
with the measured bunch profile in Fig. [3] For speed, the 
smooth lattice approximation was used. 



i = 1, 2) due to a strong coupling between the horizontal 
and vertical planes, as discussed in Sec. IIII E[ As in the 
tracking simulation, we take only collisions, intrabeam 
scattering and radiation damping into account. In total 
we have six dynamical variables: e X yi, en, Ni, where eu 



is the longitudinal emittance. 
system is given by: 



The time evolution of the 




sim. turns/10 

FIG. 7. The amount of macro particles lost (integrated over 
1000 simulation turns) due to IBS and collisions from bunch 
1 (Blue) during the simulation of the test store 8908. In total 
5 x 10 5 macro particles were simulated over 7 x 10 4 turns, 
corresponding to a bunch population of 1.056 x 10 9 and 5.26 h 
in the real machine. 



de x 



c xyi 



^xyi 



xyt 



dt Tms,xy(Ni,e X yi,eu) T la d,xy T co u(Ni, e xy i, en) 
deu eu e u 

dt TiBS,i(Ni,e xy i,eu) T ra d, 2 
dN t Ni Ni 

dt TiQs,N(Ni,e xy i,en) Tc(Ni,e xy i,en) 

Here Tibs, xy, Tibs,; are the instantaneous IBS rise 
times, and Tibs,n,Tc the instantaneous lifetimes from 
debunching and collisions. Furthermore, T co n is the emit- 
tance rise time caused by core depletion in the colli- 
sions [i^]. If the transverse distribution is gaussian, the 
collision probability is much higher for the particles in the 
core of the beam. When these particles are removed, the 
remaining ones therefore have a larger transverse emit- 
tance. Core depletion is included automatically in the 
tracking through the calculation of individual collision 
probabilities for every particle but has to be added ex- 
plicitly in the ODE system. The core depletion effect is 
discussed in detail in Ref. [42j]. 

By adding additional terms to the equations, e.g. as 
was done for beam-gas scattering in Ref. @ , the method 
can easily be expanded to account for other effects lead- 
ing to particle loss or emittance growth. 

To solve the system numerically or, in certain special 
cases, analytically we use Mathematica (43|. Our imple- 
mentation allows the kinetics of unequal beam intensities 
and emittances to be treated although we shall restrict 
ourselves to equal beams in this paper. When using the 
method of pre-calculated and interpolated values of Tibs,; 
and T\BS,xy, as explained below, the gain in execution 
speed when simulating the time evolution in a store is 
more than a factor 1000 compared with the tracking. A 
10 h store takes on the order of a second to solve on a 
normal desktop computer. 

For convenience, we used MAD-X to calculate Tibs,; and 
Tms,xy The theoretical framework 1441 is an extended 
version of the Conte-Martini model 133], which includes 
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also vertical dispersion, gaussian beam distributions are 
assumed in all planes. The influence of the IBS model on 
the simulation result is discussed in Appendix [B] where 
we make a comparison between MAD-X and the modified 
Piwinski model. The evaluation of Tibs.z and T\BS,xy 
for just one set of values Ni(t), e xy i(t), eu(t) requires a 
significant amount of computer time, so it is impractical 
to do this at runtime. 

We therefore use pre-calculated values of the rise times 
evaluated on a grid of points in the relevant region of 
(exipQ)- Since Tbs ^ , it is only necessary to 
evaluate the rise times for one typical value of JV< and 
then scale. The result is a smooth surface, which can 
be interpolated with cubic polynomials in the two emit- 
tances. This initial step, which has to be done for each 
optical configuration and beam energy under considera- 
tion, takes less than an hour of computer time on a nor- 
mal desktop machine. In the Mathematica framework, 
Tms,xy and Tibs,; are replaced by rapidly evaluated in- 
terpolating functions with no loss of precision. This is the 
key step that allows a fast interactive analysis of many 
cases of interest. Since we assume round beams, we use 
T ms,xy = ( t ibL + T ms, y )/ 2 ( above transition in a flat 
accelerator lattice the uncoupled IBS calculation gives a 
large positive growth rate in the horizontal plane and a 
small damping rate in the vertical). The calculated rise 
times from MAD-x for the initial distributions are shown 
in Table El 

The debunching effect is in principle similar to Tou- 
schek scattering, with the difference that in the standard 
formulas for the Touschek effect [45| all scattering events 
leading to a loss are assumed to occur with large momen- 
tum transfers in the center of the bucket. The debunch- 
ing losses are however diffusive in nature, since most of 
the lost particles were already very close to the separa- 
trix before the scattering. To model this we use a method 
developed by Piwinski [40|, as follows. 

If we assume that the beam distributions change 
slowly, we can approximate them as being in a steady 
state during the time step in the numerical integration. 
The instantaneous lifetime arising when the longitudi- 
nal aperture cuts the tails of the gaussian distribution is 
then, in analogy with Ref. [46j , given by 



1 6 Z 
= 2122 ex p 

Tibs, at ^TiBS,l(Ni,e xy i,eu)a s 



2*1 



(11) 



Here <5 ma x is the maximum allowed fractional energy de- 
viation AE/Eq (Eq is the reference energy) and 



7rrj c E 



(12) 



is the standard deviation of AE/Eq, with fis being the 
angular synchrotron frequency and r\ c the frequency slip 
factor. Eq. (fT2j) is derived for small oscillations, which 
limits the accuracy of the model. It should be noted 
that a similar effect exists in the transverse planes due to 
the cutoff of the physical and dynamic aperture. In the 



ideal case considered here, without significant impact of 
resonances and non-linearities, this effect is negligible. 

We determine the instantaneous partial lifetime due 
to the collisions, Tc, from the total number of particles 
removed from each beam per time: 



1 1 dNi _ Ca 

T~ c ~ ~A~^T ~ A~ 



(13) 



Here £ is given by Eq. © and a by Table Hill The reduc- 
tion factor is calculated from Eq. §5§ assuming gaussian 
distributions in the longitudinal plane. 

Finally, the emittance rise time T co u for beam i due to 
core depletion has been calculated to be 42] 



1 



/CxyiNj / re v^IP -Rtot O 



Tcoii 4V2nf3* y (e xyi + e xy] f /2 ' 
Here i = 1, 2 and j — 2 for i = 1 and vice versa. 



(14) 



V. COMPARISON BETWEEN SIMULATIONS 
AND MEASUREMENTS IN RHIC 

As an example of the detailed time evolution in a 
store, Fig. [5] shows the measured and simulated lumi- 
nosity and bunch populations for store 8908 (parameters 
are given in Table IVIj) . This store was above average in 
luminosity performance but not exceptionally good. Re- 
sults from both simulation methods (tracking and ODEs) 
are shown. A very good agreement is obtained with the 
tracking method. Using ODEs, the agreement with data 
is significantly worse. This is clearly because of the lim- 
itations in the debunching model and the assumption of 
gaussian longitudinal profiles. In the remainder of this 
section, we therefore focus on the tracking method only. 



TABLE VI. Starting parameters for the example store called 
8908. The transverse starting emittance was not measured, 
but inferred from the measured luminosity, bunch populations 
and hourglass factor using Eq. ()A9Jl . It was assumed to be 
equal for both beams and planes. 



Parameter 


Unit 


value 


Ni (Blue) 


10 a particles 


1.056 


N 2 (Yellow) 


10 9 particles 


1.004 


£ 


lO^cm-V 1 


3.0 


transv. geom. rms emittance 


mm mrad 


2.34 


FWHM of bunch length 


ns 


2.6 



In order to have more statistics, we simulated, us- 
ing the tracking method, 139 stores without stochastic 
cooling with varying bunch populations and luminosity 
from RHIC Run-7. We use the fits to the measurements, 
Eqs. (JXJ) and for comparisons as they are rather accu- 
rate during the first 3 hours. We used the same current, 
shown in Fig. 01 for all stores, which introduces an error, 
but allows us to better test the predictive ability of the 
code. 
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FIG. 8. The luminosity (top), Blue bunch intensity (middle) 
and Yellow bunch intensity (bottom) from measurements and 
two simulation methods for store 8908 in RHIC. 



To measure the goodness of the tracking simulation, 
we use two benchmark parameters, £ and ip, defined as 
the average instantaneous ratio or ratio of the integrals 
between a simulated and measured quantity U over the 
first three hours: 



1 

3h 



t=3h 



t=0h 



, t=3h u 
Jt=0h Usl 



Usimjt) 

,(*) dt 



at 



//=o>fit(t)di 



(15) 



Here U is one of the quantities C,N\,N2. Some his- 
tograms of £ and i/j f° r 139 stores are shown in Fig. |9] 
We have deselected 9 stores, where the automatically cal- 
culated fit parameters were unrealistic or the beam cur- 
rent or luminosity dropped too rapidly to possibly be 
accounted for by collisions and IBS. The mean values 
and standard deviations of £ and ip are shown in Ta- 
ble IVII1 An excellent agreement is found for the bunch 
population, while the integrated luminosity is on average 
overestimated by 13%. Given the realities of practical 
machine operation, we still consider this as a good agree- 
ment. 




FIG. 9. Comparison of tracking simulation and measurement 
of integrated luminosity (top), instantaneous Blue (middle) 
and Yellow (bottom) bunch intensity for 139 stores in RHIC 
without stochastic cooling. The parameters £ and ip are de- 
fined by Eq. (fT5|). 



TABLE VII. The mean (£) and the standard deviation crj of 
the ratio between simulated and measured quantity over 139 
stores in Run-7 without stochastic cooling in RHIC. 



Quantity 


<0 








L 


1.16 


0.056 


1.13 


0.044 


Ni 


1.02 


0.032 


1.02 


0.029 


N 2 


1.01 


0.033 


1.01 


0.031 



Several possible sources of the discrepancies between 
measured and simulated luminosity exist. Apart from 
the uncertainty resulting from the approximations in the 
simulation model, and the use of the same longitudinal 
initial conditions in all stores, variations of the machine 
optics, in particular /3* between different stores could 
give rise to errors. The logged data show that in some 
stores, the luminosity varies between the IPs at STAR 
and PHENIX, although they have nominally the same 
8* 

The interaction cross sections have uncertainties of the 
order of a few 10%. Decreasing a to 180 barn in the 
simulation changes the luminosity by 4% after 3 h. The 
cross section for mutual EMD, which was used to convert 
event rate to measured luminosity, has an uncertainty of 
not less than about 5% 1151. The uncertainties in cross 
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sections could therefore explain a significant part of the 
ovcrcstimation of the luminosity in the simulation. 

Calculating IBS rise times in the smooth lattice ap- 
proximation, which we consider to be a less accurate 
model that in this case overestimates the influence of 
IBS, improves agreement with the measured luminosity 
(if) = 1.10) while keeping the excellent agreement with 
the measured intensities (£ s» 1). This need not mean 
that the used IBS rise times calculated with the more 
detailed models are too low. It could also be that some 
other process, not included in the simulation, causes ad- 
ditional beam losses. Such unaccounted processes include 
orbit variations, non-linearities, RF noise, ground mo- 
tion, triplet quadrupole errors, current noise and the ef- 
fect of dynamic aperture. These processes can not easily 
be included in the simulation, since their strengths are 
not well known. 



VI. PREDICTIONS FOR THE LHC 

As the overall agreement with data from RHIC is very 
good, we use the tracking code to make predictions for 
the LHC, both at injection and collision energy. Run 
parameters are presented in Table[IJ There arc important 
differences with respect to RHIC. The LHC uses a single 
RF system and bunches are expected to have a gaussian 
distribution in the longitudinal plane, meaning that the 
ODE method also has the potential to work well. 

In Fig. [Til] we show simulations of the time evolution at 
the injection plateau, where bunches are kept circulating 
without colliding while the ring is filled. This is expected 
to take about 10 minutes per ring but, to cover cases 
where there may be some delay in starting the ramp, we 
have simulated 1 h. IBS rise times are longer than in 
RHIC, especially in the transverse plane (see Table ITvT) . 
The transverse emittance is expected to be 2% larger af- 
ter 20 minutes and 7% larger after 1 h. The simulated 
bunch current is 1.4% lower after 20 min and 5.2% af- 
ter 1 h because of debunching losses. The profile stays 
approximately gaussian. 

Simulations of the time evolution at collision energy 
using ODEs were done in Ref. Q. Here we redo the 
calculations, including also core depletion, and compare 
with tracking for 1-3 active IPs (the three experiments 
ALICE, ATLAS and CMS are expected to study heavy- 
ion collisions) . The results are presented in Fig. [TT] At 
collision energy the dynamics in the LHC changes signif- 
icantly, since radiation damping (see Table [Vj) becomes 
important and compensates the emittance blowup from 
IBS. As particles are removed from the beams by colli- 
sions, IBS becomes weaker while damping stays constant. 
Therefore, the emittance shrinks during the store. The 
agreement between the simulation methods is excellent. 
Losses from debunching are predicted to be negligible in 
the conditions simulated for the LHC. Beam losses from 
collisions are by far the dominant effect. 
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FIG. 10. The simulated time evolution in the LHC of the 
bunch intensity and transverse rms emittance during 1 h at 
the injection plateau without collisions. The bottom plot 
shows the longitudinal bunch profile at different times as sim- 
ulated with tracking. 



VII. CONCLUSIONS AND OUTLOOK 

We have presented measurements of the time evolution 
of the luminosity and bunch intensities in both beams 
during RHIC Run- 7 and compared with two simulation 
methods: tracking and ODEs. Tracking simulations of 
139 stores without stochastic cooling show that the pa- 
rameter £, defined as the average ratio between a simu- 
lated and measured quantity, is close to 1 with a stan- 
dard deviation of only 3% for the bunch populations. An 
average over-estimate of the measured integrated lumi- 
nosity by around 13% was found. Possible sources of the 
discrepancies include the use of the same longitudinal 
starting profile in all stores, variations in /?* , the cross 
section for 2-neutron electromagnetic dissociation (used 
to convert the measured luminosity) and the neglection 
of certain factors including RF noise and dynamic aper- 
ture. 

The ODE simulations show significant discrepancies 
with measurements in RHIC, which come from the gaus- 
sian approximation of the longitudinal bunch profile and 
the limitations in the debunching model. Including non- 
gaussian profiles in this method is not possible in the 
present framework. 

We have also made predictions of the time evolution 
in the LHC, where gaussian bunches are expected and 
the two simulation methods are in very good agreement. 
In this case, the ODE method increases the execution 
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FIG. 11. The simulated time evolution in the LHC of the lu- 
minosity, bunch intensity and transverse rms emittance dur- 
ing a 10 h store at top energy with colliding beams. Results 
are shown from the tracking method (dotted lines) and the 
ODE method (solid lines) for the cases of one (black lines), 
two (blue lines) or three (red lines) active IPs taking collisions. 
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FIG. 12. Schematic view of the collision between two bunches 
at an IP, showing the three coordinate systems used: the s- 
axis is fixed in the laboratory frame and with s — at the 
IP, while the Z\ and Z2 axes move with the origin fixed in the 
center of the two bunches. The bunch movement is indicated 
by the red arrows. The transverse coordinates of all systems 
coincide, since zero crossing angle is assumed. All distances 
are measured in the laboratory frame. 
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Appendix A: Collision probability 

To derive the phase averaged collision probability P\ in 
Eq. (j4|) we consider the movement of a particle in bunch 1 
through bunch 2 at an IP, using the coordinate systems 
defined in Fig. [T2] and consider first the case with zero 
crossing angle. The two bunches move with opposite ve- 
locities v so their centers have s = ±vt. Both centers are 
at s = at the IP at time r = 0. We write the density 
of bunch i, normalized to one, as pi{x,y,Zi). The total 
number of reactions C ac per interaction cross section a 
during a single bunch crossing is [47| : 



speed by a factor 1000. The dynamics in the LHC is 
however significantly different, since radiation damping is 
a stronger effect than IBS. This gives rise to a shrinking 
emittance at collision energy (although this benefit will 
vanish rapidly if the LHC is operated at less than full 
energy). Therefore, debunching is not an issue in the 
LHC and the collisions are the main source of particle 
losses. 

Both simulation methods are suitable as testing tools 
to optimize the luminosity. If stochastic cooling were 
to be included, the tracking could be used to make pre- 
cise predictions of the improvement in luminosity due to 
stochastic cooling in one or both beams in RHIC. 



C sc = MN 1 N 2 x 

J pi(x, y,s — VT)p 2 (x,y, s + vt) dxdydsdr, (Al) 

where M is a kinematic factor defined as 

M = - «2) 2 - Wi x « 2 ) 2 /c 2 - (A2) 

With zero crossing angle and bunch velocities |ui| = 
| #2 1 = v, we have M — 2v. 

To obtain Pi , we define a distribution function p p i for 
a single particle in bunch 1 using the Dirac ^-function. 
Since zero transverse magnetic field is assumed at the IP, 
a particle in bunch 1 with spatial coordinates (xi, yi, Zi) 
in the bunch and transverse angles (x' x , y[ ) follows a 
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straight line so 



With a similar phase averaging over <ft y , Eq. (IA6|) gives 



p pl (x,y,s- vt) = 

6(s - vt - zi)6(x - [xi + x' 1 s])S(y - [yi + y[s}). (A3) 

In the most general collision routine in the tracking 
code, p 2 is determined by sorting the particles in bunch 2 
in discrete bins along the directions (x,y,z 2 ). We as- 
sume that the transverse distributions at s — r = are 
independent, writing them as p 2x {x)p 2y {u)i an d that the 
longitudinal density p 2z is independent of x and y. The 
transverse binnings are performed at s = r = 0. 

However, we must take into account that the transverse 
distributions change along s with the optical function 
Pxy{s), which we assume to be equal in both planes in 
the proximity of the IP. In the drift section around the 
IP f} X y{s) = Plyi^ + s 2 1 Ply 2 )i where f3* y is the presumed 
minimum at the IP. At a given s close to the IP, the 
transverse distributions are wider by a factor that we 
call k(s), defined as 



k(s) 



IPxy(s) 



P. 



■'•y 




(A4) 



Inserting p v \ in Eq. (|A1[) . multiplying by a, and in- 
tegrating, gives Pi for a particle with given coordinates 
(x 1 ,x' 1 ,y 1 ,y[,z 1 ): 

* / ai+x^s s * f yi+y[s \ 

Pi = 2aN 2 I — — P2z(2s - zi) ds. 



k{s) 



«(«) 



(A5) 

In the code, the integral in Eq. (| A5|) is replaced by a sum 
over all bins that a specific particle passes through. 

This method is however slow. A faster code is obtained 
by assuming gaussian distributions in x and y. Eq. (|A5|) 
then becomes 



Pi = 2aN 2 x 



x^2y 



p 2z (2s - zi) As. 

(A6) 



Here the standard deviation of coordinate u in bunch 2 
at s = is given by yj(5* y e2u- We then change to action- 
angle variables (J x ,4>x) using Eq. ([3]) and average Pi over 
X . The xi-dependent exponential in Eq. (|A6[) becomes 



-2tt 



exp 



Jx 

C-2x 



Pxy + /3J y 



exp 



2e 



^a: \ T I Jx 

Jo 



2.r 



2e 2a 



2tt 



(A7) 



Pi = 0-^2 



exp i 



"2eo 



27r^* y ye2^e2^ 

2y02z(2s - Zx) 



1 



ds. (A8) 



Eq. (|A8p is the exact phase-averaged collision probability 
for a particle in bunch 1 with coordinates (J x , J y ,zi). We 
note that the integral represents the hour glass factor for 
a particle at a given zi. 

Furthermore we approximate the hourglass factor to 
be equal for all particles, by averaging over all zi-values. 
Changing integration variable from s to z 2 for ease of 
implementation, the global hourglass factor Rh is 



Rh — 



Pz\{zi)p z2 {z 2 ) 

(zi+z 2 ) 2 
4/3* a 



1 



dzidz 2 . 



(A9) 



It can be shown that, if p zu are gaussian, Eq. 
equivalent to the well-known formulas in Ref. [4f 
With this approximation, Pi becomes 



Pi ^aN 2 - 



Rh- 



ZnPZ yy/£2xt2y 

(A10) 

Simulations show that the errors introduced for sin- 
gle particles by using a global Rh are insignificant when 
considering the whole bunch. For the cases considered in 
this article, only negligible changes in the simulated time 
evolution of the luminosity and bunch intensity appear 
if Eq. (lAlOp is used instead of Eq. (IA5|) . while the exe- 
cution time drops by a factor 15. Therefore, we only use 
the faster method for the results in this article. 

If the bunches collide with a small crossing angle 20, we 
replace Rh by a more general luminosity reduction factor 
i?t t, which accounts for both the crossing angle and the 
hourglass effect. Such a factor is derived in Ref. [49[ for 
equal gaussian beams by integration of Eq. (|A1[) . Our 
calculation is completely analogous, with the generaliza- 
tion that we assume unequal transverse beam sizes and 
unknown longitudinal distributions. The integrations 
over (y,r) can be carried out directly yielding Eq. ([6]), 
with i?tot given by Eq. ((SJ). 



Appendix B: Comparison of IBS models 

The simulation results are relatively insensitive to the 
IBS model. As an illustration, we have calculated the 
time evolution of the quantities of interest for the LHC 
collision scheme using the modified Piwinski method for 
IBS in the ODEs, instead of mad-x. Three experiments 
were assumed active. 

In Fig.[T31 we show the ratio U m /U p , where U is the lu- 
minosity, bunch population or transverse emittance. The 
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FIG. 13. The ratio of the quantities calculated using mad-x 
for IBS to the ones where the modified Piwinski model was 
used. The time dependence was solved with the ODE method 
assuming the LHC collision scheme (numerical parameters are 
shown in Table [J, with three experiments active. 
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subscript m indicates that MAD-X was used for IBS calcu- 
lations, while p stands for the modified Piwinski model. 
The emittance shows the largest deviations (around 1% 
after 10 h), which is small compared to the overall error. 
The IBS rise times as a function of transverse emittance 
are shown for both models in Fig. [T4l 



FIG. 14. Horizontal and longitudinal IBS rise times for the 
LHC collision scheme as a function of the transverse emittance 
as calculated with mad-x and the modified Piwinski model 
(mod. Piw.). All other parameters were assumed constant 
with numerical values shown in Table [I] 



[1] H. Halm, E. Forsyth, H. Foelsche, M. Harrison, 
J. Kewisch, G. Parzen, S. Peggs, E. Raka, A. Ruggiero, 
A. Stevens, S. Tepikian, P. Thieberger, D. Trbojevic, 
J. Wei, E. Willen, S. Ozaki, and S. Y. Lee, Nucl. Instrum. 
Methods Phys. Res., Sect. A 499, 245 (2003). 

[2] O. S. Briining, P. Collier, P. Lebrun, S. Myers, R. Ostojic, 
J. Poole, and P. Proudlock (editors), CERN-2004-003- 
VI (2004). 

[3] K. Hubner and E. Keil, IEEE Transactions on Nuclear 

Science NS-32, 1632 (1985). 
[4] D. Brandt, K. Eggert, and A. Morsch, CERN SL/94-04 

(AP)(1994). 

[5] A. J. Baltz, M. J. Rhoades-Brown, and J. Weneser, Phys. 
Rev. E 54, 4233 (1996). 

[6] J. M. Jowett, H. H. Braun, M. I. Gresham, E. Man- 
ner, A. N. Nicholson, and E. Shaposhnikova, Proc. of 
the European Particle Accelerator Conf. 2004, Lucerne, 
578(2004). 

[7] J. Wei and A. Ruggiero, BNL-45269 AD/RHIC-81(1990). 
[8] K. Gounder, J. Marriner, and S. Mishra, Proc. of the 

Particle Accelerator Conf. 2003, Portland, 3434(2003). 
[9] A. O. Sidorin, I. N. Meshkov, I. A. Seleznev, 



A. V. Smirnov, E. M. Syresin, and G. V. Trubnikov, 

INucl. Inst, fc Methods A 558, 325 (2006). 
[10] S. Y. Zhang, RHIC C-A/AP/297(2007). 
[11] A. Drees et al, Proc. of the Particle Accelerator Conf. 

2007, Albuquerque, New Mexico, USA, 722(2007). 
[12] M. Blaskiewicz and J. M. Brennan, Proceedings of COOL 

2007, Bad Kreuznach, Germany, 125(2007). 
[13] M. Blaskiewicz and J. M. Brennan, Phys. Rev. ST Accel. 

Beams 10, 061001 (Jun 2007). 
[14] M. Blaskiewicz, J. M. Brennan, and F. Severino, Phys. 

Rev. Letters 100, 174802 (2008). 
[15] A. J. Baltz, C. Chasman, and S. N. White, Nucl. Instrum. 

Methods Phys. Res., Sect. A 417, 1 (1998). 
[16] C. Adler, A. Denisov, E. Garcia, 

M. Murray, H. Stroebele, and S . White, 

|Nucl. Instrum. Methods Phys. Res., Sect.~A] 470, 

488 (2001), ISSN 0168-9002. 
[17] D. Trbojevic, RHIC/AP/136(1997). 

[18] D. Trbojevic, H. C. Hsueh, W. MacKay, A. Drees, and 
R. Fliller, in Proceedings of the 2001 Particle Accelerator 
Conference, Chicago (2001) p. 3141. 

[19] J. M. Jowett, Proceedings of the LHC Performance 



15 



Workshop, Chamonix XIII, 55(2005). 
[20] J. M. Jowett, Proceedings of the LHC Performance [35 

Workshop, Chamonix XIV, 88(2005). 
[21] A.W. Chao, M. Tigner (editors), Handbook of Accelerator 

Physics and Engineering (World Scientific, 1998). [36 
[22] |http : //cern . ch/mad7| 
[23] S. R. Klein, Nucl. Inst. & Methods A 459, 51 (2001). [37 
[24] J. M. Jowett, J. B. Jeanneret, and K. Schindl, Proc. of the [38 

Particle Accelerator Conf. 2003, Portland, 1682(2003). 
[25] R. Bruce, D. Bocian, S. Gilardoni, and J. M. Jowett, [39 

|Phys. Rev. ST Accel. Beams| 12, 071002 (Jul 2009). 
[26] H. Meier, Z. Halabuka, K. Hencken, D. Trautmann, and [40 

G. Baur, Phys. Rev. A 63, 032713 (2001). 
[27] I. A. Pshenichnov, J. P. Bondorf, I . N. Mishustin, A. Ven- [41 

tura, and S. Masetti, |Phys. Rev. C| 64, 024903 (Jul 2001). 
[28] ALICE Collaboration: F. Carminati, P. Foka, [42 

P. Giubellino, A. Morsch, G. Paic, J.-P. Revol, [43 

K. Safarik, Y. Schutz, and U. A. Wiedemann (editors), [44 

Journal of Physics G: Nuclear and Particle Physics 30, [45 

1517 (2004). " [46 
[29] A. J. Baltz, BNL-66644(1999). 
[30] J. Dunlop, private communication. [47 
[31] A. Piwinski, Proc. of the 9th International Conference on 

High Energy Accelerators, Stanford, CA, 405(1974). [ 
[32] J. Bjorken and S. Mtingwa, Particle Accelerators 13, 115 

(1983). [49 
[33] M. Conte and M. Martini, Part. Acc. 17, 1 (1985). 
[34] K. Kubo and K. Oide, Phys. Rev. ST Accel. Beams 4, 



124401 (2001). 

K. L. F. Bane, H. Hayano, K. Kubo, T. Naito, T. Okugi, 
and J. Urakawa, Phys. Rev. ST Accel. Beams 5, 084403 
(2002). 

S. Nagaitsev, Phys. Rev. ST Accel. Beams 8, 064403 
(2005). 

F. Zimmermann, private communication. 

R. Siemann, AIP conference 127, AIP New York, 
368(1985). 

G. Parzen, Proc. of the European Particle Accelerator 
Conf. 1988, Rome, Italy, 821(1988). 

M. Blaskiewicz, Proc. of the Particle Accelerator Conf. 

2003, Kno xville, Tennessee, USA, 310(2005) 

C. Kim, in Proceedings of the 1997 Particle Accelerator Conference , 

Vol. 1 (1997) pp . 790-792 vol.1. 

R. Br uce(2009), ar Xiv:0911.5627vll 

http : //www .wolf ram. com 

F. Zimmermann, CERN-AB-2006-002(2005). 

A. Piwinski, DESY 98-179(1998). 

A. Piwinski, Proc. of the CERN Accelerator School, Gif- 
sur-Yvette, France CERN-85-19 (1985). 

C. M0ller, K. Danske Vidensk. Selsk. Mat.-Fys. Medd. 
23 (1945). 

M. A. Furman, Proc. of the Particle Accelerator Conf. 
1991, San Fransisco, California, 422(1991). 

B. Muratori, CERN LHC Project Note 301(2002). 



